*Purpose: Post-Matching Estimation with Difference-in-Difference
** By state regressions for manuscript tables, including binary and continuous treatment ** 
** Also by PA regressions for main manuscript tables, including binary and continuous treatment **

clear


*******************************************
** Run from Master file for consistency ***
*******************************************

*ssc install schemepack, replace
net install grc1leg,from( http://www.stata.com/users/vwiggins/)

**0. Load data
use "$data\MASTER_MATCHED_R2.dta", clear
	rename cell_uniqueID cellID
    xtset cellID post
	
**1. Define globals 
global invarcontrols elevation_km slope_100deg NearCity_dist_100km NearCity_population i.LandSurface_code i.Bioclimate_code i.Geology_code
global timevarcontrols Precipitation Temperature


** 2. Estimations by state

	
if 1{
* Binary treatment	
	xtreg NetOverNF post treat_post $timevarcontrols , fe i(cellID) cluster(cd_mun)
	est store all, title(All)
	
	foreach v in Bahia EspiritoSanto MinasGerais Parana RioDeJaneiro SaoPaulo {
		display "***`v'***"
		xtreg NetOverNF post treat_post $timevarcontrols if state=="`v'", fe i(cellID) cluster(cd_mun)
		est store `v', title(`v')
	}
	
	*Print table in stata
	esttab all Bahia EspiritoSanto Parana RioDeJaneiro SaoPaulo MinasGerais, keep(treat_post post) order(treat_post post) mtitles("Full" "Bahia" "Esp Santo" "Parana" "Rio de Jan" "Sao Paulo" "Min Gerais")
	
	*LaTeX Output Tables
	esttab all Bahia EspiritoSanto Parana RioDeJaneiro SaoPaulo MinasGerais using "$tables\TableByStateBinary.tex", replace label starlevels(* 0.10 ** 0.05 *** 0.01)  se nonotes ///
	keep(post treat_post ) ///
	order(treat_post post ) ///
	varwidth(3)  style(tex) mtitles("Full" "Bahia" "Esp Santo" "Parana" "Rio de Jan" "Sao Paulo" "Min Gerais") b(%12.3f) se(%12.3f)
	
	*Coefplot
	coefplot (all, label(All) ciopts(lwidth(0.75) lcolor("51 166 166")) msymbol(Oh) msize(large) mcolor("51 166 166")) ///
	(Bahia, label(Bahia) ciopts(lwidth(0.75) lcolor("242 149 68")) msymbol(O) msize(large) mcolor("242 149 68")) ///
	(Parana, label(Parana) ciopts(lwidth(0.75) lcolor("242 80 65")) msymbol(Sh) msize(large) mcolor("242 80 65")) ///
	(EspiritoSanto, label(Espirito Santo) ciopts(lwidth(0.75) lcolor("50 54 115")) msymbol(S) msize(large) mcolor("50 54 115")) ///
	(RioDeJaneiro, label(Rio De Janeiro) ciopts(lwidth(0.75) lcolor("140 31 102")) msymbol(Dh) msize(large) mcolor("140 31 102")) ///
	(SaoPaulo, label(Sao Paulo) ciopts(lwidth(0.75) lcolor("166 94 51")) msymbol(D) msize(large) mcolor("166 94 51"))  ///
	(MinasGerais, label(Minas Gerais) ciopts(lwidth(0.75) lcolor("41 51 51")) msymbol(Th) msize(large) mcolor("41 51 51")),  ///
	vertical yline(0) keep(treat_post) coeflabels(treat_post="Treatment Effect") legend(pos(6) rows(1)) xlabel(, noticks nolabel) scheme(white_tableau) ytitle("Marginal effect on forest restoration") ylabel(-0.1(0.1)0.3 ) 

	graph export "${figures}\Figure3.eps", as(eps) replace
	*graph export "${figures}\Figure3_1000mSpillover.eps", as(eps) replace
	*graph export "${figures}\Figure3_NoPAs.eps", as(eps) replace
	
	*Coefplot
	coefplot (Bahia, label(Bahia) ciopts(lwidth(0.75) lcolor("51 166 166")) msymbol(O) msize(large) mcolor("51 166 166")) ///
	(Parana, label(Paraná) ciopts(lwidth(0.75) lcolor("242 80 65")) msymbol(Sh) msize(large) mcolor("242 80 65")) ///
	(EspiritoSanto, label(Espírito Santo) ciopts(lwidth(0.75) lcolor("50 54 115")) msymbol(S) msize(large) mcolor("50 54 115")) ///
	(RioDeJaneiro, label(Rio de Janeiro) ciopts(lwidth(0.75) lcolor("140 31 102")) msymbol(Dh) msize(large) mcolor("140 31 102")) ///
	(SaoPaulo, label(São Paulo) ciopts(lwidth(0.75) lcolor("166 94 51")) msymbol(D) msize(large) mcolor("166 94 51"))  ///
	(MinasGerais, label(Minas Gerais) ciopts(lwidth(0.75) lcolor("41 51 51")) msymbol(Th) msize(large) mcolor("41 51 51")), ///
	vertical yline(0) keep(treat_post) coeflabels(treat_post="Treatment Effect") legend(pos(12) rows(1)) xlabel( , noticks nolabel)  scheme(white_tableau) ytitle("Marginal effect on forest restoration") saving("${figures}\ByState", replace)
	
	gr combine "${figures}\ByState" "${figures}\Fines2", col(1) imargin(0 0 0 0) 
	gr export "${figures}\Figure6.eps", replace	
	
	gr combine "${figures}\ByState" "${figures}\Fines", col(1) imargin(0 0 0 0) 
	gr export "${figures}\FigureCoeffFinesCapita.eps", replace

	
estimates clear
	
* Continuous treatment

quietly: xtreg NetOverNF post post_treat_percent $timevarcontrols , fe i(cellID) cluster(cd_mun)
	est store all, title(All)
	
foreach v in Bahia EspiritoSanto MinasGerais Parana RioDeJaneiro SaoPaulo {
	quietly: xtreg NetOverNF post post_treat_percent $timevarcontrols if state=="`v'", fe i(cellID) cluster(cd_mun)

	est store `v', title(`v')
}

	*LaTeX Output Tables
	esttab all Bahia EspiritoSanto Parana RioDeJaneiro SaoPaulo MinasGerais using "$tables\TableByStateContinuous.tex", replace label starlevels(* 0.10 ** 0.05 *** 0.01)  se nonotes ///
	keep(post post_treat_percent ) ///
	order(post_treat_percent post ) ///
	varwidth(3)  style(tex) mtitles("Full" "Bahia" "Esp Santo" "Parana" "Rio de Jan" "Sao Paulo" "Min Gerais") b(%12.3f) se(%12.3f)
	
estimates clear 
}

if 0{

** 3. Estimations with PA

foreach v in Bahia EspiritoSanto MinasGerais Parana RioDeJaneiro SaoPaulo {
	quietly: xtreg NetOverNF post treat_post PA_post PA_treatpost $timevarcontrols  if state=="`v'", fe i(cellID) cluster(cd_mun)

	est store `v', title(`v')
}
xtreg NetOverNF post treat_post PA_post PA_treatpost $timevarcontrols, fe i(cellID) cluster(cd_mun)
	est store Full, title(Full)
	
	*LaTeX Output Tables
	esttab  Full Bahia EspiritoSanto Parana RioDeJaneiro SaoPaulo MinasGerais using "$tables\PAbyState.tex", replace label starlevels(* 0.10 ** 0.05 *** 0.01)  se nonotes ///
	keep(post treat_post PA_post PA_treatpost  ) ///
	order(treat_post post PA_post PA_treatpost  ) ///
	varwidth(3)  style(tex) mtitles("Full" "Bahia" "Esp Santo" "Parana" "Rio de Jan" "Sao Paulo" "Min Gerais") b(%12.3f) se(%12.3f)
		
	graph export "${figures}\Figure4.png", replace
}

** 4. Analyzing heterogeneity associated with our theoretical model 

if 0 {
	g proximity = 1/NearCity_dist_100km

	
	xtreg NetOverNF i.post##c.treat_percent##c.proximity $timevarcontrols , fe i(cellID) cluster(cd_mun)
	margins, dydx(treat_percent) at(post=(1) proximity=( 10.89 4.45 3.21 2.16 1.18))  post noestimcheck
	
	coefplot, ciopts(lwidth(0.75) lcolor("242 149 68")) msymbol(O) msize(large) mcolor("242 80 65") ///
	coeflabels(1._at treat_percent = "9 km" 2._at treat_percent = "23 km" 3._at treat_percent = "31 km" 4._at treat_percent = "46 km" 5._at treat_percent = "83 km") ///
	xtitle("Marginal effect of restoration, proportional treatment")
		graph export "${figures}\FigureDist.eps", as(eps) replace
		

graph bar (median) NearCity_dist_100km, label over(state) ytitle("Km to nearest city") 
graph export "${figures}\MedianDist.png", replace

	

}